#### NOTE: WHEN FINISHED REPLACE soclib WITH cultlib ?
####       ALSO SEARCH FOR "SOCIAL" AND REPLACE AS NEEDE

#library(rio)
library(car)
library(rdrobust)
library(AER)
library(corrplot)
library(sensemakr)

## reading in raw data and formatting
load("CulturalLiberalismDataJOP.RData")

# this function standardizes each item (also making it a vector)
scale.nomatrix <- function(x){as.vector(scale(x))}

## Creating data frame of traditional values indicators
## note: signs flipped for some vars so higher values
##       always correspond to more traditional values
## Cultural Liberalism is the reverse of traditional values
DAT$itv1 <- -scale.nomatrix(DAT$q1_6)
DAT$itv2 <- -scale.nomatrix(DAT$q1_7)
DAT$itv3 <- -scale.nomatrix(DAT$q1_8)
DAT$itv4 <- -scale.nomatrix(DAT$q1_9)
DAT$itv5 <- -scale.nomatrix(DAT$q1_10)
DAT$itv6 <- scale.nomatrix(DAT$q1_11)
DAT$itv7 <- scale.nomatrix(DAT$q1_12)
DAT$itv8 <- scale.nomatrix(DAT$q6_1)
DAT$itv9 <- -scale.nomatrix(DAT$q8_1)
DAT$itv10 <- -scale.nomatrix(DAT$q8_2)
DAT$itv11 <- -scale.nomatrix(DAT$q8_3)
DAT$itv12 <- scale.nomatrix(DAT$q8_4)
DAT$itv13 <- -scale.nomatrix(DAT$q9)
DAT$itv14 <- scale.nomatrix(DAT$q10)
DAT$itv15 <- -scale.nomatrix(DAT$q11)
DAT$itv16 <- scale.nomatrix(DAT$neigh_drug)
DAT$itv17 <- scale.nomatrix(DAT$neigh_race)
DAT$itv18 <- scale.nomatrix(DAT$neigh_AIDS)
DAT$itv19 <- scale.nomatrix(DAT$neigh_immig)
DAT$itv20 <- scale.nomatrix(DAT$neigh_homo)
DAT$itv21 <- scale.nomatrix(DAT$neigh_relig)
DAT$itv22 <- scale.nomatrix(DAT$neigh_couples)
DAT$itv23 <- scale.nomatrix(DAT$neigh_language)
DAT$itv24 <- -scale.nomatrix(DAT$q13_1)
DAT$itv25 <- -scale.nomatrix(DAT$q13_2)
DAT$itv26 <- -scale.nomatrix(DAT$q13_3)
DAT$itv27 <- -scale.nomatrix(DAT$q13_4)
DAT$itv28 <- -scale.nomatrix(DAT$q13_5)
DAT$itv29 <- -scale.nomatrix(DAT$q13_6)

## data frame of (raw, with NAs) indicators for trad. values
TVAL <- cbind.data.frame(DAT$itv1, DAT$itv2, DAT$itv3, DAT$itv4, 
                         DAT$itv5, DAT$itv6, DAT$itv7, DAT$itv8,
                         DAT$itv9, DAT$itv10, DAT$itv11, DAT$itv12,
                         DAT$itv13, DAT$itv14, DAT$itv15,
                         DAT$itv16, DAT$itv17, DAT$itv18, DAT$itv19,
                         DAT$itv20, DAT$itv21, DAT$itv22, DAT$itv23,
                         DAT$itv24, DAT$itv25, DAT$itv26, DAT$itv27,
                         DAT$itv28, DAT$itv29)

# esitmating missing-imputed PCA
library(missMDA)
# imputing  missing values
tradval.imputePCA <- imputePCA(TVAL)
# estimating PCA with imputed data
library(FactoMineR)
pca.tradval <- PCA(tradval.imputePCA$completeObs)

# creating traditional values scale
# based on PCA first component with imputation
# Note: making it negative b/c it's flipped
#       so in DAT, higher values = more traditional
DAT$tradval <- pca.tradval$ind$coord[,1]
DAT$tradval <- (DAT$tradval - mean(DAT$tradval)) / sd(DAT$tradval)
# creating social liberalism variable
# that's just negative of tradval
DAT$soclib <- -DAT$tradval

# recoding childhood SES recall question
DAT$childhoodSES <- recode(DAT$p22, "'Lower class'=1; 'Working class'=2; 'Lower middle class'=3; 'Upper middle class'=4; 'Upper class'=5; else=NA")
DAT$fatheredu <- DAT$p20

# creating full dataset
DAT.WITH2019 <- DAT
DAT.ONLY2019 <- subset(DAT.WITH2019, GRADYEAR==2019)
# subsetting data to drop GRADYEAR==2019
DAT <- subset(DAT, GRADYEAR!=2019)

## assessing compliance
## table of proportion attending university
## among those above/below bac passage threshold
prop.table(table(DAT$average >= 6, DAT$p6), margin=1)


#############
## TABLE 1 ##
#############
# ITT RD analyses
# predicting pre-treatment covaraites
## pre-treatment covariate analyses
## (all ITT)
## father's education
lm.fatheredu <- lm(fatheredu ~ average>=6,
                   data=DAT)
summary(lm.fatheredu)
confint(lm.fatheredu)
##
lm.childhoodSES <- lm(childhoodSES ~ average>=6,
                      data=DAT)
summary(lm.childhoodSES)
confint(lm.childhoodSES)
## humanities track
lm.humanist <- lm(humanist ~ average>=6,
                  data=DAT)
summary(lm.humanist)
confint(lm.humanist)
## urban HS
lm.HS_urban <- lm(HS_urban=="MUNICIPIU" ~ average>=6,
                  data=DAT)
summary(lm.HS_urban)
confint(lm.HS_urban)


## within .15 bandwidth
lm.fatheredu.15 <- lm(fatheredu ~ average>=6,
                      data=subset(DAT, abs(average-6)<=.15))
summary(lm.fatheredu.15)
confint(lm.fatheredu.15)
##
lm.childhoodSES.15 <- lm(childhoodSES ~ average>=6,
                         data=subset(DAT, abs(average-6)<=.15))
summary(lm.childhoodSES.15)
confint(lm.childhoodSES.15)
## humanities track
lm.humanist.15 <- lm(humanist ~ average>=6,
                     data=subset(DAT, abs(average-6)<=.15))
summary(lm.humanist.15)
confint(lm.humanist.15)
## urban HS
lm.HS_urban.15 <- lm(HS_urban=="MUNICIPIU" ~ average>=6,
                     data=subset(DAT, abs(average-6)<=.15))
summary(lm.HS_urban.15)
confint(lm.HS_urban.15)

## within .1 bandwidth
lm.fatheredu.1 <- lm(fatheredu ~ average>=6,
                     data=subset(DAT, abs(average-6)<=.1))
summary(lm.fatheredu.1)
confint(lm.fatheredu.1)
##
lm.childhoodSES.1 <- lm(childhoodSES ~ average>=6,
                        data=subset(DAT, abs(average-6)<=.1))
summary(lm.childhoodSES.1)
confint(lm.childhoodSES.1)
## humanities track
lm.humanist.1 <- lm(humanist ~ average>=6,
                    data=subset(DAT, abs(average-6)<=.1))
summary(lm.humanist.1)
confint(lm.humanist.1)
## urban HS
lm.HS_urban.1 <- lm(HS_urban=="MUNICIPIU" ~ average>=6,
                    data=subset(DAT, abs(average-6)<=.1))
summary(lm.HS_urban.1)
confint(lm.HS_urban.1)

## within .05 bandwidth
lm.fatheredu.05 <- lm(fatheredu ~ average>=6,
                      data=subset(DAT, abs(average-6)<=.05))
summary(lm.fatheredu.05)
confint(lm.fatheredu.05)
##
lm.childhoodSES.05 <- lm(childhoodSES ~ average>=6,
                         data=subset(DAT, abs(average-6)<=.05))
summary(lm.childhoodSES.05)
confint(lm.childhoodSES.05)
## humanities track
lm.humanist.05 <- lm(humanist ~ average>=6,
                     data=subset(DAT, abs(average-6)<=.05))
summary(lm.humanist.05)
confint(lm.humanist.05)
## urban HS
lm.HS_urban.05 <- lm(HS_urban=="MUNICIPIU" ~ average>=6,
                     data=subset(DAT, abs(average-6)<=.05))
summary(lm.HS_urban.05)
confint(lm.HS_urban.05)

## closest values to the cutoff
lm.fatheredu.closest <- lm(fatheredu ~ average>=6,
                           data=subset(DAT, average==6 | average==5.983333))
summary(lm.fatheredu.closest)
confint(lm.fatheredu.closest)
## childhood SES
lm.childhoodSES.closest <- lm(childhoodSES ~ average>=6,
                              data=subset(DAT, average==6 | average==5.983333))
summary(lm.childhoodSES.closest)
confint(lm.childhoodSES.closest)
## humanities
lm.humanist.closest <- lm(humanist ~ average>=6,
                          data=subset(DAT, average==6 | average==5.983333))
summary(lm.humanist.closest)
confint(lm.humanist.closest)
## urban HS
lm.HS_urban.closest <- lm(HS_urban=="MUNICIPIU" ~ average>=6,
                          data=subset(DAT, average==6 | average==5.983333))
summary(lm.HS_urban.closest)
confint(lm.HS_urban.closest)

## CONTINUITY-BASED RD FOR PRE-TREATMENT COVARIATES
## Father's Education
clustervar <- DAT$average
rd.fatheredu <- rdrobust(y=DAT$fatheredu, x=DAT$average, 
                         c=6, h=.2,
                         vce = "hc0", cluster = clustervar)
summary(rd.fatheredu)
## Childhood SES
clustervar <- DAT$average
rd.childhoodSES <- rdrobust(y=DAT$childhoodSES, x=DAT$average, 
                            c=6, h=.2,
                            vce = "hc0", cluster = clustervar)
summary(rd.childhoodSES)
## Humanities
clustervar <- DAT$average
rd.humanist <- rdrobust(y=DAT$humanist, x=DAT$average, 
                        c=6, h=.2,
                        vce = "hc0", cluster = clustervar)
summary(rd.humanist)
## HS urban
clustervar <- DAT$average
rd.HS_urban <- rdrobust(y=DAT$HS_urban=="MUNICIPIU", x=DAT$average, 
                        c=6, h=.2,
                        vce = "hc0", cluster = clustervar)
summary(rd.HS_urban)




##############
## FIGURE 3 ##
##############
# RD plot of attending university and cultural liberalism
# by bac score (each point is a unique bac score value
# with size of points proportional to sample sizes)
par(mfrow=c(1,2))
#
plot(tapply(DAT$p6, DAT$average, mean, na.rm=TRUE) ~ sort(unique(DAT$average)),
     cex=table(DAT$average)/15, pch=20, bty="n", ylim=c(0,1),
     xlab="Bac Score", ylab="Proportion attending university")
abline(v=6, lty=3, col="red")
abline(h=tapply(DAT$p6, DAT$average>=6, mean, na.rm=TRUE), lty=3, col="gray")
#
plot(tapply(DAT$soclib, DAT$average, mean, na.rm=TRUE) ~ sort(unique(DAT$average)), 
     cex=table(DAT$average)/15,
     pch=20, bty="n", 
     xlab="Bac Score", ylab="Average Cultural Liberalism")
abline(v=6, lty=3, col="red")



#############
## TABLE 2 ##
#############
# local randomization RD estimates
# Panel A: First Stage
lm.firststage <- lm(p6 ~ average>=6, 
                    data=subset(DAT,
                                !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage)
confint(lm.firststage)
# first stage .15 
lm.firststage.15 <- lm(p6 ~ average>=6, 
                       data=subset(DAT, abs(average-6)<=.15 &
                                     !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.15)
confint(lm.firststage.15)
## first stage .1
lm.firststage.1 <- lm(p6 ~ average>=6, 
                      data=subset(DAT, abs(average-6)<=.1 &
                                    !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.1)
confint(lm.firststage.1)
## first stage .05
lm.firststage.05 <- lm(p6 ~ average>=6, 
                       data=subset(DAT, abs(average-6)<=.05 & 
                                     !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.05)
confint(lm.firststage.05)
## first stage closest 
lm.firststage.closest <- lm(p6 ~ average>=6, 
                            data=subset(DAT, (average==6 | average==5.983333) &
                                          !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.closest)
confint(lm.firststage.closest)

# Panel B: ITT
lm.itt <- lm(soclib ~ average>=6, 
             data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt)
confint(lm.itt)
## ITT results .15
lm.itt.15 <- lm(soclib ~ average>=6, 
                data=subset(DAT, abs(average-6)<=.15 & 
                              !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.15)
confint(lm.itt.15)
## ITT results .10
lm.itt.10 <- lm(soclib ~ average>=6, 
                data=subset(DAT, abs(average-6)<=.10 & 
                              !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.10)
confint(lm.itt.10)
## ITT results .05
lm.itt.05 <- lm(soclib ~ average>=6, 
                data=subset(DAT, abs(average-6)<=.05 & 
                              !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.05)
confint(lm.itt.05)
## ITT closest 
lm.itt.closest <- lm(soclib ~ average>=6, 
                     data=subset(DAT, (average==6 | average==5.983333) & 
                                   !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.closest)
confint(lm.itt.closest)

# Panel C: Effect of University Attendance
lm.tsls <- ivreg(soclib ~ p6 | average>=6,
                 data=DAT)
summary(lm.tsls)
confint(lm.tsls)
# tsls .15
lm.tsls.15 <- ivreg(soclib ~ p6 | average>=6,
                    data=subset(DAT, abs(average-6)<=.15))
summary(lm.tsls.15)
confint(lm.tsls.15)
# tsls .1
lm.tsls.1 <- ivreg(soclib ~ p6 | average>=6,
                   data=subset(DAT, abs(average-6)<=.1))
summary(lm.tsls.1)
confint(lm.tsls.1)
# tsls .05
lm.tsls.05 <- ivreg(soclib ~ p6 | average>=6,
                    data=subset(DAT, abs(average-6)<=.05))
summary(lm.tsls.05)
confint(lm.tsls.05)
# tsls closest
lm.tsls.closest <- ivreg(soclib ~ p6 | average>=6,
                         data=subset(DAT, average==6 | average==5.983333))
summary(lm.tsls.closest)
confint(lm.tsls.closest)



#############
## TABLE 3 ##
#############
# continuity-based RD analysis
# Panel A: First Stage
clustervar <- DAT$average
rd.firststage <- rdrobust(y=DAT$p6[!is.na(DAT$average) & 
                                     !is.na(DAT$p6) & 
                                     !is.na(DAT$soclib)], 
                          x=DAT$average[!is.na(DAT$average) & 
                                          !is.na(DAT$p6) & 
                                          !is.na(DAT$soclib)], 
                          c=6, h=.2,
                          vce = "hc0", 
                          cluster = clustervar[!is.na(DAT$average) & 
                                                 !is.na(DAT$p6) & 
                                                 !is.na(DAT$soclib)])
summary(rd.firststage)
# Panel B: ITT
clustervar <- DAT$average
rd.itt <- rdrobust(y=DAT$soclib[!is.na(DAT$average) & 
                                  !is.na(DAT$p6) & 
                                  !is.na(DAT$soclib)], 
                   x=DAT$average[!is.na(DAT$average) & 
                                   !is.na(DAT$p6) & 
                                   !is.na(DAT$soclib)],
                   c=6, h=.2,
                   vce = "hc0", 
                   cluster = clustervar[!is.na(DAT$average) & 
                                          !is.na(DAT$p6) & 
                                          !is.na(DAT$soclib)])
summary(rd.itt)
# Panel C: Effect of University Attendance
clustervar <- DAT$average
rd.main <- rdrobust(y=DAT$soclib, x=DAT$average, 
                    c=6, h=.2, fuzzy=DAT$p6,
                    vce = "hc0", cluster = clustervar)
summary(rd.main)



#############
## TABLE 4 ##
#############
# creating labeled matrix of tradval indicators
TVAL.LAB <- TVAL
colnames(TVAL.LAB) <- c("accept.homosexuality",
                        "accept.prostitution",
                        "accept.abortion",
                        "accept.divorce",
                        "accept.premarital",
                        "accept.beatwife",
                        "accept.beatkids",
                        "churchmember",
                        "jobs.men",
                        "jobs.native",
                        "woman.earn",
                        "woman.jobind",
                        "attend.relig",
                        "importance.god",
                        "world.better",
                        "neighbor.addict",
                        "neighbor.otherrace",
                        "neighbor.aids",
                        "neighbor.immigrant",
                        "neighbor.homosexual",
                        "neighbor.otherrelig",
                        "neighbor.unmarried",
                        "neighbor.otherlang",
                        "parents.proud",
                        "mother.work",
                        "men.betterleaders",
                        "educ.boys",
                        "men.betterexecutives",
                        "housewife.fulfill")

# creating subcategory scores
TVAL.IMPUTED <- as.data.frame(tradval.imputePCA$completeObs)
colnames(TVAL.IMPUTED) <- colnames(TVAL.LAB)
DAT.WITH2019$soclib.sub.gender <- -scale(TVAL.IMPUTED$accept.beatwife +
                                       TVAL.IMPUTED$accept.beatkids + 
                                       TVAL.IMPUTED$jobs.men + 
                                       TVAL.IMPUTED$woman.earn + 
                                       TVAL.IMPUTED$woman.jobind + 
                                       TVAL.IMPUTED$mother.work + 
                                       TVAL.IMPUTED$men.betterleaders +
                                       TVAL.IMPUTED$educ.boys +
                                       TVAL.IMPUTED$men.betterexecutives + 
                                       TVAL.IMPUTED$housewife.fulfill)
DAT$soclib.sub.gender <- DAT.WITH2019$soclib.sub.gender[DAT.WITH2019$GRADYEAR!=2019]
#
DAT.WITH2019$soclib.sub.relig <- -scale(TVAL.IMPUTED$churchmember +
                                      TVAL.IMPUTED$attend.relig +
                                      TVAL.IMPUTED$importance.god +
                                      TVAL.IMPUTED$world.better +
                                      TVAL.IMPUTED$neighbor.otherrelig)
DAT$soclib.sub.relig <- DAT.WITH2019$soclib.sub.relig[DAT.WITH2019$GRADYEAR!=2019]
#
DAT.WITH2019$soclib.sub.moral <- -scale(TVAL.IMPUTED$accept.homosexuality +
                                      TVAL.IMPUTED$accept.prostitution +
                                      TVAL.IMPUTED$accept.abortion + 
                                      TVAL.IMPUTED$accept.divorce + 
                                      TVAL.IMPUTED$accept.premarital + 
                                      TVAL.IMPUTED$neighbor.addict + 
                                      TVAL.IMPUTED$neighbor.aids + 
                                      TVAL.IMPUTED$neighbor.homosexual + 
                                      TVAL.IMPUTED$neighbor.unmarried + 
                                      TVAL.IMPUTED$parents.proud)
DAT$soclib.sub.moral <- DAT.WITH2019$soclib.sub.moral[DAT.WITH2019$GRADYEAR!=2019]
#
DAT.WITH2019$soclib.sub.race <- -scale(TVAL.IMPUTED$jobs.native +
                                     TVAL.IMPUTED$neighbor.otherrace +
                                     TVAL.IMPUTED$neighbor.immigrant +
                                     TVAL.IMPUTED$neighbor.otherlang)
DAT$soclib.sub.race <- DAT.WITH2019$soclib.sub.race[DAT.WITH2019$GRADYEAR!=2019]

## Table of correlations of traditional values subscores:
cor(cbind(DAT$soclib.sub.gender, DAT$soclib.sub.moral, DAT$soclib.sub.race, DAT$soclib.sub.relig))
# significance tests for correlations
cor.test(DAT$soclib.sub.gender, DAT$soclib.sub.moral)
cor.test(DAT$soclib.sub.gender, DAT$soclib.sub.race)
cor.test(DAT$soclib.sub.gender, DAT$soclib.sub.relig)
cor.test(DAT$soclib.sub.moral, DAT$soclib.sub.race)
cor.test(DAT$soclib.sub.moral, DAT$soclib.sub.relig)
cor.test(DAT$soclib.sub.race, DAT$soclib.sub.relig)



#############
## TABLE 5 ##
#############
# Local RD analyses of subscores
# gender
tsls.gender <- ivreg(DAT$soclib.sub.gender ~ DAT$p6 | DAT$average>=6)
summary(tsls.gender)
confint(tsls.gender)
# moral
tsls.moral <- ivreg(DAT$soclib.sub.moral ~ DAT$p6 | DAT$average>=6)
summary(tsls.moral)
confint(tsls.moral)
# race
tsls.race <- ivreg(DAT$soclib.sub.race ~ DAT$p6 | DAT$average>=6)
summary(tsls.race)
confint(tsls.race)
# relig
tsls.relig <- ivreg(DAT$soclib.sub.relig ~ DAT$p6 | DAT$average>=6)
summary(tsls.relig)
confint(tsls.relig)
#
# Continuity-based Subscore-specific RDs
clustervar <- DAT$average
rd.gender <- rdrobust(y=DAT$soclib.sub.gender, x=DAT$average, 
                      c=6, h=.2, fuzzy=DAT$p6,
                      vce = "hc0", cluster = clustervar)
summary(rd.gender)
#
clustervar <- DAT$average
rd.moral <- rdrobust(y=DAT$soclib.sub.moral, x=DAT$average, 
                     c=6, h=.2, fuzzy=DAT$p6,
                     vce = "hc0", cluster = clustervar)
summary(rd.moral)
#
clustervar <- DAT$average
rd.race <- rdrobust(y=DAT$soclib.sub.race, x=DAT$average, 
                    c=6, h=.2, fuzzy=DAT$p6,
                    vce = "hc0", cluster = clustervar)
summary(rd.race)
#
clustervar <- DAT$average
rd.relig <- rdrobust(y=DAT$soclib.sub.relig, x=DAT$average, 
                     c=6, h=.2, fuzzy=DAT$p6,
                     vce = "hc0", cluster = clustervar)
summary(rd.relig)



#############################
## APPENDIX TABLES/FIGURES ##
#############################

##############
## Table B1 ##
##############
# gender
prop.table(table(DAT.WITH2019$p11))
# single (never married)
prop.table(table(DAT.WITH2019$p12=="Single"))
# one or more children
prop.table(table(DAT.WITH2019$p13>=1))
# employed
prop.table(table(DAT.WITH2019$p14_1=="Yes"))
# Romanian ethnicity
prop.table(table(DAT.WITH2019$p23[DAT.WITH2019$p23!=""]=="român"))
# father's education level (4 cat)
prop.table(table(DAT.WITH2019$fatheredu4))
# childhood SES
prop.table(table(DAT.WITH2019$childhoodSES))
# current SES
prop.table(table(DAT.WITH2019$p15[DAT.WITH2019$p15!=""]))
# age
prop.table(table(2020-DAT.WITH2019$p17))


########################
## APPENDIX FIGURE C1 ##
########################
corrplot(cor(TVAL.LAB, use="p"), method="ellipse")


#######################
## APPENDIX TABLE C1 ##
#######################
# NOTE: variables are in different order in appendix table
cbind(names(TVAL.LAB), pca.tradval$var$coord[,1])



########################
## APPENDIX FIGURE C2 ##
########################
# density plot of cultural liberalism with overall average
# and average for failers, passers
plot(density(DAT$soclib), 
     xlim=c(-2,2), lwd=2, bty="n", main="",
     yaxt="n", xlab="cultural liberalism")
abline(v=mean(DAT$soclib[DAT$average>=6], na.rm=TRUE),
       lty=2, lwd=2)
abline(v=mean(DAT$soclib[DAT$average<6], na.rm=TRUE),
       lty=3, lwd=2)
legend('topright',
       c("average for passers", "average for failers"),
       lwd=2, lty=2:3)
text(1.4, 0, "more culturally liberal ->")
text(-1.4, 0, "<- less culturally liberal")



###############
## FIGURE D1 ##
###############
bandwidth.radius <- c(.075,.1, .125,.15, .175,.2)
bandwidth.est.table <- matrix(NA, nrow=length(bandwidth.radius), ncol=3)
rownames(bandwidth.est.table) <- bandwidth.radius
colnames(bandwidth.est.table) <- c("tau_hat", "ci_lower", "ci_upper")

for(i in 1:length(bandwidth.radius)){
  foo <- rdrobust(y=DAT$soclib,
                  x=DAT$average, 
                  c=6, h=bandwidth.radius[i], 
                  fuzzy=DAT$p6,
                  vce = "hc0", 
                  cluster = DAT$average)
  bandwidth.est.table[i, ] <- c(foo$coef[1,], foo$ci[3,])
}

plot(bandwidth.radius, bandwidth.est.table[,1], 
     xlab="Bandwidth Radius", 
     ylab="Estimated effect of university attendance",
     bty="n", ylim=c(-.5, 1.5), 
     pch=20, col="red", xaxt="n")
axis(1, at=bandwidth.radius, labels=substr(bandwidth.radius, 2, 5), las=2)
for(i in 1:length(bandwidth.radius)){
  segments(bandwidth.radius[i], bandwidth.est.table[i,2], 
           bandwidth.radius[i], bandwidth.est.table[i,3])
}
abline(h=0, col="gray", lty=3)



###############
## FIGURE D2 ##
###############
donut.radius <- c(0, .02, .04, .06, .08, .1)
donut.est.table <- matrix(NA, nrow=length(donut.radius), ncol=3)
rownames(donut.est.table) <- donut.radius
colnames(donut.est.table) <- c("tau_hat", "ci_lower", "ci_upper")

for(i in 1:length(donut.radius)){
  foo <- rdrobust(y=DAT$soclib[DAT$average < 6 - donut.radius[i] |
                                 DAT$average >= 6 + donut.radius[i]],
                  x=DAT$average[DAT$average < 6 - donut.radius[i] | 
                                  DAT$average >= 6 + donut.radius[i]], 
                  c=6, h=.2, 
                  fuzzy=DAT$p6[DAT$average < 6 - donut.radius[i] | 
                                 DAT$average >= 6 + donut.radius[i]],
                  vce = "hc0", 
                  cluster = DAT$average[DAT$average < 6 - donut.radius[i] | 
                                          DAT$average >= 6 + donut.radius[i]])
  donut.est.table[i, ] <- c(foo$coef[1,], foo$ci[3,])
}

plot(donut.radius, donut.est.table[,1], 
     xlab="Donut Hole Radius", ylab="Estimated effect of university attendance",
     bty="n", ylim=range(donut.est.table), pch=20, col="red")
for(i in 1:length(donut.radius)){
  segments(donut.radius[i], donut.est.table[i,2], 
           donut.radius[i], donut.est.table[i,3])
}
abline(h=0, col="gray", lty=3)




##############
## TABLE D1 ##
##############
# local randomization estimates
# by discipline
# Humanist:
# Panel A: First Stage
lm.firststage.humanist <- lm(p6 ~ average>=6, 
                             data=subset(DAT, humanist==1 & 
                                           !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.humanist)
confint(lm.firststage.humanist)
# Panel B: ITT
lm.itt.humanist <- lm(soclib ~ average>=6, 
                      data=subset(DAT, humanist==1 & 
                                    !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.humanist)
confint(lm.itt.humanist)
# Panel C: Effect of University Attendance
lm.tsls.humanist <- ivreg(soclib ~ p6 | average>=6,
                          data=subset(DAT, humanist==1))
summary(lm.tsls.humanist)
confint(lm.tsls.humanist)
#
## SCIENCE
# Panel A: First Stage
lm.firststage.science <- lm(p6 ~ average>=6, 
                            data=subset(DAT, humanist==0 & 
                                          !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.science)
confint(lm.firststage.science)
# Panel B: ITT
lm.itt.science <- lm(soclib ~ average>=6, 
                     data=subset(DAT, humanist==0 & 
                                   !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.science)
confint(lm.itt.science)
# Panel C: Effect of University Attendance
lm.tsls.science <- ivreg(soclib ~ p6 | average>=6,
                         data=subset(DAT, humanist==0))
summary(lm.tsls.science)
confint(lm.tsls.science)


##############
## TABLE D2 ##
##############
# Continuity-based RD estimates by discipline
# Note: using only observations with non-missing
# values on all variables so will match main analysis
#
# Humanist:
# Panel A: First Stage
clustervar <- DAT$average
rd.firststage.humanist <- rdrobust(y=DAT$p6[DAT$humanist==1 & 
                                              !is.na(DAT$average) & 
                                              !is.na(DAT$p6) & 
                                              !is.na(DAT$soclib)], 
                                   x=DAT$average[DAT$humanist==1 & 
                                                   !is.na(DAT$average) & 
                                                   !is.na(DAT$p6) & 
                                                   !is.na(DAT$soclib)], 
                                   c=6, h=.2,
                                   vce = "hc0", 
                                   cluster = clustervar[DAT$humanist==1 & 
                                                          !is.na(DAT$average) & 
                                                          !is.na(DAT$p6) & 
                                                          !is.na(DAT$soclib)])
summary(rd.firststage.humanist)
# Panel B: ITT
clustervar <- DAT$average
rd.itt.humanist <- rdrobust(y=DAT$soclib[DAT$humanist==1 & 
                                           !is.na(DAT$average) & 
                                           !is.na(DAT$p6) & 
                                           !is.na(DAT$soclib)], 
                            x=DAT$average[DAT$humanist==1 & 
                                            !is.na(DAT$average) & 
                                            !is.na(DAT$p6) & 
                                            !is.na(DAT$soclib)],
                            c=6, h=.2,
                            vce = "hc0", 
                            cluster = clustervar[DAT$humanist==1 & 
                                                   !is.na(DAT$average) & 
                                                   !is.na(DAT$p6) & 
                                                   !is.na(DAT$soclib)])
summary(rd.itt.humanist)
# Panel C: Effect of University Attendance
clustervar <- DAT$average
rd.humanist <- rdrobust(y=DAT$soclib[DAT$humanist==1], 
                        x=DAT$average[DAT$humanist==1], 
                        c=6, h=.2, fuzzy=DAT$p6[DAT$humanist==1],
                        vce = "hc0", cluster = clustervar[DAT$humanist==1])
summary(rd.humanist)
#
# Sciences:
# Panel A: First Stage
clustervar <- DAT$average
rd.firststage.nonhumanist <- rdrobust(y=DAT$p6[DAT$humanist==0 & 
                                                 !is.na(DAT$average) & 
                                                 !is.na(DAT$p6) & 
                                                 !is.na(DAT$soclib)], 
                                      x=DAT$average[DAT$humanist==0 & 
                                                      !is.na(DAT$average) & 
                                                      !is.na(DAT$p6) & 
                                                      !is.na(DAT$soclib)], 
                                      c=6, h=.2,
                                      vce = "hc0", 
                                      cluster = clustervar[DAT$humanist==0 & 
                                                             !is.na(DAT$average) & 
                                                             !is.na(DAT$p6) & 
                                                             !is.na(DAT$soclib)])
summary(rd.firststage.nonhumanist)
# Panel B: ITT
clustervar <- DAT$average
rd.itt.nonhumanist <- rdrobust(y=DAT$soclib[DAT$humanist==0 & 
                                              !is.na(DAT$average) & 
                                              !is.na(DAT$p6) & 
                                              !is.na(DAT$soclib)], 
                               x=DAT$average[DAT$humanist==0 & 
                                               !is.na(DAT$average) & 
                                               !is.na(DAT$p6) & 
                                               !is.na(DAT$soclib)],
                               c=6, h=.2,
                               vce = "hc0", 
                               cluster = clustervar[DAT$humanist==0 & 
                                                      !is.na(DAT$average) & 
                                                      !is.na(DAT$p6) & 
                                                      !is.na(DAT$soclib)])
summary(rd.itt.nonhumanist)
# Panel C: Effect of University Attendance
clustervar <- DAT$average
rd.nonhumanist <- rdrobust(y=DAT$soclib[DAT$humanist==0], 
                           x=DAT$average[DAT$humanist==0], 
                           c=6, h=.2, fuzzy=DAT$p6[DAT$humanist==0],
                           vce = "hc0", cluster = clustervar[DAT$humanist==0])
summary(rd.nonhumanist)


##############
## TABLE D3 ##
##############
# Covariate adjusted RD
#
# Local Randomization RD wtih Covariates
# Panel A: First Stage
lm.firststage.covs <- lm(p6 ~ I(average>=6) + childhoodSES + 
                           fatheredu + humanist + 
                           I(HS_urban == "MUNICIPIU"), 
                         data=subset(DAT,
                                     !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.covs)
confint(lm.firststage.covs)
# Panel B: ITT
lm.itt.covs <- lm(soclib ~ I(average>=6) + childhoodSES + 
                    fatheredu + humanist + 
                    I(HS_urban == "MUNICIPIU"), 
                  data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.covs)
confint(lm.itt.covs)
# Panel C: Effect of University Attendance
iv.covs <- ivreg(soclib ~ p6 + childhoodSES + 
                   fatheredu + humanist + 
                   I(HS_urban == "MUNICIPIU") | 
                   I(average>=6) + childhoodSES + 
                   fatheredu + humanist + 
                   I(HS_urban == "MUNICIPIU"), 
                 data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(iv.covs)
confint(iv.covs)
#
# Continuity-Based RD with Covariates
#
clustervar <- DAT$average
Z <- cbind(DAT$fatheredu, DAT$childhoodSES,
           DAT$humanist, DAT$HS_urban=="MUNICIPIU")
colnames(Z) <- c("fatheredu", "childhoodSES",
                 "humanist", "HS_urban")
covs.missing <- is.na(DAT$fatheredu) | is.na(DAT$childhoodSES) |
  is.na(DAT$humanist) | is.na(DAT$HS_urban) |
  is.na(DAT$p6) | is.na(DAT$soclib) | is.na(DAT$average)
# first stage with covs
rd.firststage.covs <- rdrobust(y=DAT$p6[!covs.missing], 
                               x=DAT$average[!covs.missing], 
                               covs=Z[!covs.missing, ],
                               c=6, h=.2,
                               vce = "hc0", 
                               cluster = clustervar[!covs.missing])
summary(rd.firststage.covs)
# ITT with covs
rd.itt.covs <- rdrobust(y=DAT$soclib[!covs.missing], 
                        x=DAT$average[!covs.missing], 
                        covs=Z[!covs.missing, ],
                        c=6, h=.2,
                        vce = "hc0", 
                        cluster = clustervar[!covs.missing])
summary(rd.itt.covs)
# Main Estimate with covs
rd.main.covs <- rdrobust(y=DAT$soclib[!covs.missing], 
                         x=DAT$average[!covs.missing], 
                         covs=Z[!covs.missing, ],
                         c=6, h=.2, fuzzy=DAT$p6[!covs.missing],
                         vce = "hc0", cluster = clustervar[!covs.missing])
summary(rd.main.covs)



###############
## FIGURE D3 ##
###############
## Effect plot for all Cultural values Indicators
CVAL.INDICATORS <- -cbind.data.frame(DAT$itv1, DAT$itv2, DAT$itv3, DAT$itv4, 
                                     DAT$itv5, DAT$itv6, DAT$itv7, DAT$itv8,
                                     DAT$itv9, DAT$itv10, DAT$itv11, DAT$itv12,
                                     DAT$itv13, DAT$itv14, DAT$itv15,
                                     DAT$itv16, DAT$itv17, DAT$itv18, DAT$itv19,
                                     DAT$itv20, DAT$itv21, DAT$itv22, DAT$itv23,
                                     DAT$itv24, DAT$itv25, DAT$itv26, DAT$itv27,
                                     DAT$itv28, DAT$itv29)
# making all 0 to 1
for(i in 1:ncol(CVAL.INDICATORS)){ 
  CVAL.INDICATORS[,i] <- (CVAL.INDICATORS[,i] - min(CVAL.INDICATORS[,i], na.rm=TRUE)) / 
    max(CVAL.INDICATORS[,i] - min(CVAL.INDICATORS[,i], na.rm=TRUE), na.rm=TRUE)
}
# cultural values indocator names
colnames(CVAL.INDICATORS) <- c("homosexuality", "prostitution",
                               "abortion", "divorce", 
                               "premarital sex", "beat wife",
                               "beat children", 
                               "church member", 
                               "men jobs", "native jobs",
                               "women earn", "woman job",
                               "church attend", "god important",
                               "science better", 
                               "neighbor drugs",
                               "neighbor other race",
                               "neighbor AIDS",
                               "neighbor immigrant",
                               "neighbor gay",
                               "neighbor other religion",
                               "neighbor unmarried",
                               "neighbor other language",
                               "make parents proud",
                               "mother work",
                               "men better leaders",
                               "univ. boys",
                               "men better exec.",
                               "housewife fulfilling"
)
# estimating continuity-based RD for each indicator
CVAL.EFFECTS <- matrix(NA, nrow=ncol(CVAL.INDICATORS), ncol=3)
colnames(CVAL.EFFECTS) <- c("estimate", "2.5%", "97.5%")
library(rdrobust)
for(i in 1:nrow(CVAL.EFFECTS)){
  foo.rd <- rdrobust(y=CVAL.INDICATORS[,i], x=DAT$average, 
                     c=6, h=.2, fuzzy=DAT$p6,
                     vce = "hc0", cluster = DAT$average)
  CVAL.EFFECTS[i, ] <- c(foo.rd$coef[1], foo.rd$ci[3,])
}
# plot of cultural values indicator effects
par(mar=c(2,2,2,4))
plot(c(NA, NA, NA, CVAL.EFFECTS[1:29, 1]),
     bty="n", xlab="", xaxt="n", xlim=c(4, 33),
     ylab="", ylim=c(-.5,.5))
axis(side= 4, at=seq(-.4, .4, by=.2))
mtext("Estimated Effect", side=4, line=2)
for(i in 1:nrow(CVAL.EFFECTS[1:29,])){
  segments(i+3, CVAL.EFFECTS[1:29,][i, 2], i+3, CVAL.EFFECTS[1:29,][i,3], lwd=2)
  text(i+3, -.5, colnames(CVAL.INDICATORS[1:29])[i], 
       adj=c(0, .5), srt=90, cex=.8)
}
abline(h=0, lty=3)
points(c(NA, NA, NA, CVAL.EFFECTS[1:29, 1]), pch=19)



###############
## FIGURE D4 ##
###############
# Sensitivity Analysis Contour Plots
# for Models Without Covariates
# re-running first stage model
lm.firststage <- lm(p6 ~ average>=6, 
                    data=subset(DAT,
                                !is.na(average) & !is.na(p6) & !is.na(soclib)))
# re-running reduced form model
lm.itt <- lm(soclib ~ average>=6, 
             data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
# sensitivity analysis of first stage
firststage.sens <- sensemakr(model = lm.firststage, 
                             treatment = "average >= 6TRUE",
                             kd = 2)
reducedform.sens <- sensemakr(model = lm.itt, 
                              treatment = "average >= 6TRUE",
                              kd = 2)
#
par(mfrow=c(1,2))
#
plot(reducedform.sens, 
     sensitivity.of = "t-value", 
     lim = 0.2, lim.y = 0.2, nlevels = 6,
     xlab="Partyal R2 of confounder(s) with the instrument",
     main="Reduced Form Sensitivity (no covariates)")
#
plot(firststage.sens, 
     sensitivity.of = "t-value", 
     lim = 0.6, lim.y = 0.6, nlevels = 6,
     xlab="Partial R2 of confounder(s) with the instrument",
     ylab="Partial R2 of confounder(s) with the treatment",
     main="First Stage Sensitivity (no covariates)")




################
## FIGURE D5 ##
################
# Sensitivity Analysis Contour Plots
# for Models With Covariates
# first stage model with covariates
lm.firststage.covs <- lm(p6 ~ I(average>=6) + childhoodSES + 
                           fatheredu + humanist + 
                           I(HS_urban == "MUNICIPIU"), 
                         data=subset(DAT,
                                     !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.firststage.covs)
confint(lm.firststage.covs)
# reduced form model with covariates
lm.itt.covs <- lm(soclib ~ I(average>=6) + childhoodSES + 
                    fatheredu + humanist + 
                    I(HS_urban == "MUNICIPIU"), 
                  data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(lm.itt.covs)
confint(lm.itt.covs)
# IV regression with covs
iv.covs <- ivreg(soclib ~ p6 + childhoodSES + 
                   fatheredu + humanist + 
                   I(HS_urban == "MUNICIPIU") | 
                   I(average>=6) + childhoodSES + 
                   fatheredu + humanist + 
                   I(HS_urban == "MUNICIPIU"), 
                 data=subset(DAT, !is.na(average) & !is.na(p6) & !is.na(soclib)))
summary(iv.covs)
confint(iv.covs)
# sensitivity analysis of first stage
firststage.sens.covs <- sensemakr(model = lm.firststage.covs, 
                                  treatment = "I(average >= 6)TRUE",
                                  benchmark_covariates=names(coef(lm.itt.covs))[-(1:2)],
                                  kd = 10)
reducedform.sens.covs <- sensemakr(model = lm.itt.covs, 
                                   treatment = "I(average >= 6)TRUE",
                                   benchmark_covariates=names(coef(lm.itt.covs))[-(1:2)],
                                   kd = 10)
#
par(mfrow=c(1,2))
#
plot(reducedform.sens.covs, 
     sensitivity.of = "t-value", 
     lim = 0.2, lim.y = 0.2, nlevels = 6,
     xlab="Partial R2 of confounder(s) with the instrument",
     main="Reduced Form Sensitivity (with covariates)")
#
plot(firststage.sens.covs, 
     sensitivity.of = "t-value", 
     lim = 0.6, lim.y = 0.6, nlevels = 6,
     xlab="Partial R2 of confounder(s) with the instrument",
     ylab="Partial R2 of confounder(s) with the treatment",
     main="First Stage Sensitivity (with covariates)")



###############
## FIGURE E1 ##
###############
year.est.table <- matrix(NA, nrow=5, ncol=3)
rownames(year.est.table) <- 2015:2019
colnames(year.est.table) <- c("tau_hat", "ci_lower", "ci_upper")
year.locrand.table <- year.est.table
# looping over graduation years in our sample (including 2019)
for(i in 1:5){
  subyear <- (2015:2019)[i]
  # conducting same RD analysis as our main result
  # using only gradyear-specific data
  foo <- rdrobust(y=DAT.WITH2019$soclib[DAT.WITH2019$GRADYEAR==subyear],
                  x=DAT.WITH2019$average[DAT.WITH2019$GRADYEAR==subyear], 
                  c=6, h=.2, fuzzy=DAT.WITH2019$p6[DAT.WITH2019$GRADYEAR==subyear],
                  vce = "hc0", cluster = DAT.WITH2019$average[DAT.WITH2019$GRADYEAR==subyear])
  # storing estimated coef and CI
  year.est.table[i, ] <- c(foo$coef[1,], foo$ci[3,])
  # note: using "standard" estimate and "robust" CI
  foo.locrand <- ivreg(DAT.WITH2019$soclib[DAT.WITH2019$GRADYEAR==subyear]  ~ 
                         DAT.WITH2019$p6[DAT.WITH2019$GRADYEAR==subyear] | 
                         DAT.WITH2019$average[DAT.WITH2019$GRADYEAR==subyear] >= 6)
  year.locrand.table[i, ] <- c(foo.locrand$coefficients[2], 
                               confint(foo.locrand)[2,1],
                               confint(foo.locrand)[2,2])
}
# creating GRADYEAR sub-analyses plot
plot((2015:2019)-.1, year.est.table[,1], 
     xlab="Year of Graduation", xaxt="n",
     ylab="Estimated effect of university attendance",
     bty="n", xlim=c(2014.8, 2019.2), ylim=range(year.est.table), 
     pch=20, col="red")
points((2015:2019)+.1, year.locrand.table[,1],
       pch=20, col="blue")
axis(1, at=2015:2019, labels=2015:2019)
for(i in 1:5){
  segments(((2015:2019)-.1)[i], year.est.table[i,2], ((2015:2019)-.1)[i], year.est.table[i,3])
  segments(((2015:2019)+.1)[i], year.locrand.table[i,2], ((2015:2019)+.1)[i], year.locrand.table[i,3])
}
abline(h=0, col="gray", lty=3)
legend('bottomleft', c("continutiy-based", "local randomization"),
       pch=20, col=c("red", "blue"))



##############
## TABLE I1 ##
##############
# RD Estimates of Effect of University Attendance
# by High School Location
# 
rd.nonurban <- rdrobust(y=DAT$soclib[DAT$HS_urban!="MUNICIPIU"], 
                        x=DAT$average[DAT$HS_urban!="MUNICIPIU"], 
                        c=6, h=.2,
                        fuzzy=DAT$p6[DAT$HS_urban!="MUNICIPIU"],
                        vce = "hc0", 
                        cluster = clustervar[DAT$HS_urban!="MUNICIPIU"])
summary(rd.nonurban)
#
clustervar <- DAT$average
rd.urban <- rdrobust(y=DAT$soclib[DAT$HS_urban=="MUNICIPIU"], 
                     x=DAT$average[DAT$HS_urban=="MUNICIPIU"], 
                     c=6, h=.2,
                     fuzzy=DAT$p6[DAT$HS_urban=="MUNICIPIU"],
                     vce = "hc0", 
                     cluster = clustervar[DAT$HS_urban=="MUNICIPIU"])
summary(rd.urban)



##############
## TABLE I2 ##
##############
# RD Estimates of Effect of University Attendance
# by Academic Concentration
clustervar <- DAT$average
rd.science <- rdrobust(y=DAT$soclib[DAT$humanist==0], 
                       x=DAT$average[DAT$humanist==0], 
                       c=6, h=.2,
                       fuzzy=DAT$p6[DAT$humanist==0],
                       vce = "hc0", 
                       cluster = clustervar[DAT$humanist==0])
summary(rd.science)
# 
rd.humanities <- rdrobust(y=DAT$soclib[DAT$humanist==1], 
                          x=DAT$average[DAT$humanist==1], 
                          c=6, h=.2,
                          fuzzy=DAT$p6[DAT$humanist==1],
                          vce = "hc0", 
                          cluster = clustervar[DAT$humanist==1])
summary(rd.humanities)



##############
## TABLE I3 ##
##############
# Continuity-based Subscore-specific RDs
clustervar <- DAT$average
rd.gender <- rdrobust(y=DAT$soclib.sub.gender, x=DAT$average, 
                      c=6, h=.2, fuzzy=DAT$p6,
                      vce = "hc0", cluster = clustervar)
summary(rd.gender)
#
clustervar <- DAT$average
rd.moral <- rdrobust(y=DAT$soclib.sub.moral, x=DAT$average, 
                     c=6, h=.2, fuzzy=DAT$p6,
                     vce = "hc0", cluster = clustervar)
summary(rd.moral)
#
clustervar <- DAT$average
rd.race <- rdrobust(y=DAT$soclib.sub.race, x=DAT$average, 
                    c=6, h=.2, fuzzy=DAT$p6,
                    vce = "hc0", cluster = clustervar)
summary(rd.race)
#
clustervar <- DAT$average
rd.relig <- rdrobust(y=DAT$soclib.sub.relig, x=DAT$average, 
                     c=6, h=.2, fuzzy=DAT$p6,
                     vce = "hc0", cluster = clustervar)
summary(rd.relig)
# Local RD analyses of subscores
# gender
tsls.gender <- ivreg(DAT$soclib.sub.gender ~ DAT$p6 | DAT$average>=6)
summary(tsls.gender)
confint(tsls.gender)
# gender .15
tsls.gender.15 <- ivreg(soclib.sub.gender ~ p6 | average>=6,
                        data=subset(DAT, abs(average-6)<=.15))
summary(tsls.gender.15)
confint(tsls.gender.15)
# gender .10
tsls.gender.10 <- ivreg(soclib.sub.gender ~ p6 | average>=6,
                        data=subset(DAT, abs(average-6)<=.10))
summary(tsls.gender.10)
confint(tsls.gender.10)
# gender .05
tsls.gender.05 <- ivreg(soclib.sub.gender ~ p6 | average>=6,
                        data=subset(DAT, abs(average-6)<=.05))
summary(tsls.gender.05)
confint(tsls.gender.05)
# gender closest
tsls.gender.closest <- ivreg(soclib.sub.gender ~ p6 | average>=6,
                             data=subset(DAT, average==6 | average==5.983333))
summary(tsls.gender.closest)
confint(tsls.gender.closest)
# moral
tsls.moral <- ivreg(DAT$soclib.sub.moral ~ DAT$p6 | DAT$average>=6)
summary(tsls.moral)
confint(tsls.moral)
# moral .15
tsls.moral.15 <- ivreg(soclib.sub.moral ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.15))
summary(tsls.moral.15)
confint(tsls.moral.15)
# moral .10
tsls.moral.10 <- ivreg(soclib.sub.moral ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.10))
summary(tsls.moral.10)
confint(tsls.moral.10)
# moral .05
tsls.moral.05 <- ivreg(soclib.sub.moral ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.05))
summary(tsls.moral.05)
confint(tsls.moral.05)
# moral closest
tsls.moral.closest <- ivreg(soclib.sub.moral ~ p6 | average>=6,
                            data=subset(DAT, average==6 | average==5.983333))
summary(tsls.moral.closest)
confint(tsls.moral.closest)
# race
tsls.race <- ivreg(DAT$soclib.sub.race ~ DAT$p6 | DAT$average>=6)
summary(tsls.race)
confint(tsls.race)
# race .15
tsls.race.15 <- ivreg(soclib.sub.race ~ p6 | average>=6,
                      data=subset(DAT, abs(average-6)<=.15))
summary(tsls.race.15)
confint(tsls.race.15)
# race .10
tsls.race.10 <- ivreg(soclib.sub.race ~ p6 | average>=6,
                      data=subset(DAT, abs(average-6)<=.10))
summary(tsls.race.10)
confint(tsls.race.10)
# race .05
tsls.race.05 <- ivreg(soclib.sub.race ~ p6 | average>=6,
                      data=subset(DAT, abs(average-6)<=.05))
summary(tsls.race.05)
confint(tsls.race.05)
# race closest
tsls.race.closest <- ivreg(soclib.sub.race ~ p6 | average>=6,
                           data=subset(DAT, average==6 | average==5.983333))
summary(tsls.race.closest)
confint(tsls.race.closest)
# relig
tsls.relig <- ivreg(DAT$soclib.sub.relig ~ DAT$p6 | DAT$average>=6)
summary(tsls.relig)
confint(tsls.relig)
# relig .15
tsls.relig.15 <- ivreg(soclib.sub.relig ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.15))
summary(tsls.relig.15)
confint(tsls.relig.15)
# relig .10
tsls.relig.10 <- ivreg(soclib.sub.relig ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.10))
summary(tsls.relig.10)
confint(tsls.relig.10)
# relig .05
tsls.relig.05 <- ivreg(soclib.sub.relig ~ p6 | average>=6,
                       data=subset(DAT, abs(average-6)<=.05))
summary(tsls.relig.05)
confint(tsls.relig.05)
# relig closest
tsls.relig.closest <- ivreg(soclib.sub.relig ~ p6 | average>=6,
                            data=subset(DAT, average==6 | average==5.983333))
summary(tsls.relig.closest)
confint(tsls.relig.closest)
